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Abstract 

In this work, the relationship between splines and the control theory has been an- 
alyzed. We show that spline functions can be constructed naturaly from the control 
theory. By establishing a framework based on control theory, we provide a simple 
and systematic way to construct splines. We have constructed the traditional spline 
functions including the polynomial splines and the classical exponential spline. We 
have also discovered some new spline functions such as trigonometric splines and the 
combination of polynomial, exponential and trigonometric splines. The method pro- 
posed in this paper is easy to implement. Some numerical experiments are performed 
to investigate properties of different spline approximations. 


1. Introduction. 

Spline functions are well known and are widely used for practical approximation of func- 
tions or more commonly for fitting smooth curves through preassigned points. Spline tech- 
niques have the advantage over most approximation and interpolation techniques in that 
they are computational feasible. Most of the published spline algorithms are for polynomial 
splines and the vast preponderance are for cubic splines. There is a small but excellent lit- 
erature on the so called exponential splines and there is an even smaller literature on splines 
with more or less arbitrary nodal functions, [9, 3]. 

In this paper we will present a common frame work for splines that includes polynomial 
splines of all orders and generalized exponential splines of all orders. This common frame 
work is based on ideas from linear control theory. Let’s recall some basic ideas from control 
theory. A linear control system is a differential equation 

= Ax(t) + Bu(t) 
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where x G R n , u G R m and the matrices A and B are constant matrices of compatible 
dimension. The vector x is the state of the system and the vector u is the control. The idea 
is that we can use the control u to steer the state from point to point in the state space R n . 

We can think of the first component of x as representing the position of the system and for 
appropriate A the second coordinate is the velocity, the third acceleration, etc. A common 
situation, for example in air traffic control, is to specific the position that the system must 
be in at a sequence of times. So in fact what we have is a set of points through which 
the system must traverse at specified times. One could fit these points with a spline curve 
and then ask for the control that would move the system along that trajectory. In fact this 
can be done but we will show that the control law can be developed from natural control 
theoretic principles that will move the system through the points at the desired times and the 
resulting curve will be piecewise analytic and will have 2n — 1 continuous derivatives, i.e. a 
generalized spline. With this framework we can construct a wide variety of spline functions. 

If the matrix A is nilpotent then the resulting construction is just that for polynomial splines. 

If the matrix is 2 x 2 and one eigenvalue is zero and the other is a nonzero real number then 
the spline is the usual exponential spline. In general the nodal functions are the coordinate 
functions of the matrix function e At . 

In this paper we give a unified treatment of all of the common one dimensional spline 
functions using simple ideas from control theory. It is coming to be understood that there is 
a large overlap between linear control theory and elementary numerical analysis. Eigenvalue 
methods are know to be closely related to the theory of the matrix Riccati equation [2], there 
are close relations between observability and quadrature techniques [8], system identification 
and Prony’s method are very similar [1] and now we see that the spline constructions and 
basic linear controllability are manifestations of the same phenomena. 

In Section 2 we review basic material from the theory of linear control systems that is 
needed for the development and give a condition that characterizes the optimal control law 
that generates the spline functions. In Section 3 we give the details of the construction 
of spline functions using control theory and in Section 4 we classify the possible classes of 
spline functions that arise from the control theoretic construction. Iir Section 5 we examine % 
in detail some of the particular classes from Section 4 and finally in Section 6 we present a m 
series of numerical examples comparing the various classes. 
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2. Some results from the control theory. 


In this section we collect a series of results from linear control theory. Most can be found 
in any control theory textbook. See, for example, the book by Brockett, [5]. 

Consider the linear system: 

d 


dt 


x(t) = Ax{t) + £>u(*), t € [0, T], 
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and the observation function 

y(t) = <fx(t ), £ = (1,0, •■•,0). 

Let us divide [0, T] into n subintervals as 

o = to < h < • • • < t„_ 1 <t n = T, 


( 2 . 1 ) 


( 2 . 2 ) 


(2.3) 


and define hk = tk — h-i, the length of the rth subinterval. Our goal is to find a control 
law u 6 C m-2 [0,T] that drives the system (2.1) from x(0) = x° to x(T ) = x T such that the 
observation function y(t) satisfies the interpolation conditions 

y{tk) = ctk, k= 1, •••,« — 1. ^ (2.4) 

Furthermore, u(t) minimizes the functional 

f u(s) 2 ds. (2-5) 

Jo 

Such a control is called an optimal control. 

Definition. The system (2.1) is called controllable if for any x ° , x 1 ' , and r > 0, there is 
a u(f) such that, 

r =~x(t) ± +* [ T e A{T ~ s) bu(3)ds. 

Jo 

Theorem 2.1 : The system (2.1) is controllable if and only if 

rank (b, Ab, • • • , A m-1 6) = m. (2-6) 
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For the special matrix A as in (2.2), it is easy to verify that 


(b,Ab,.--,A n 



0 

0 

... o 

1 ^ 

0 
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... 1 

• • 
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1 

... * 

* 

U 

* 

... * 

* / 


(2.7) 


and hence the condition (2.6) is satisfied. From Theorem 2.1, the system (2.1) is controllable. 

Theorem 2.2 : The system (2.1) is controllable, if and only if the matrix J e A3 bb r e A 3 ds 
is invertible. 


For the special matrix in (2.2), we then define 


M(t) = ([ e- A3 bb r e- ATs ds)~ 1 . (2.8) 

Theorem 2.3 : When the system (2.1) is controllable, a control that moves the system from 
x(t) = p L to x(t) = p R given by 

u(t) = P'e~ ATt (Jt e~ Aa bb r e~ AT3 ds)~ l {e~ M p R - e~ At -p L ), (2.9) 

minimizes the functional J(y ) = J v 2 (s)ds among all controls that move the system from 
x(t) = p L to x(t) = pR. 

Theorem 2.4 : When the system (2.1) is controllable, a control u (E C m ~ 2 [0,T] that moves 
the system from «(0) = x°, passing through (? x(t k ) = oik, to x(T) = xT is given by 

n — 1 m 

«(<) = E &/*(<) + E7.-tf(0. (2 - 10) 

*= i >=i 


with 


f m _ r t < t k , 

/jk(t) - { o t>tk, k = 1,. 


,n - 1, 


9i {t) = <? e A ^-% i = !,• 


where ef = (1, 0, • • • , 0), • • • , e£ = (0, • • • , 0, 1), and 0 k ’s, Ti ’ are determined byn-1 interpo- 
lation conditions ^xftk) = a k and m boundary conditions x(T) = x? . Moreover, the control 
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(2.10) minimizes the functional J(v) = j v 2 (s)ds among all functions v £ C m 2 [0 ,T] that 
drives the system (2.1) from f(0) = x°, passing through <?x{tk) = a k , to x(T) = X 1 * . 

The construction of the optimal control is based on Hilbert space techniques and is based 
on writing the interior constraints in terms of a linear variety defined in terms of the functions 
fk{t) and the terminal constraints in terms of the functions gi(t ). Once the constraints are 
written in terms of linear varieties the form of the optimal control is clear based on the 
orthogonal complement of the intersection of the varieties, this is a standard technique and 
is found for example in [7] or [6]. The proof of Theorem 2.4 is based on two facts: (i) u(t) 
defined by (2.10) is m — 2 times continuously differentiable; (ii) /j t’s and gfs are n + m — 1 
linearly independent functions. In the following, we verify (i) and (ii). 

(i) From the construction (2.10), we only need to show that 

fk\tk) = 0> r = 0, • • • , m 2, k = !,••*, n — 1. 

Indeed, for all k = 1, • • • , n — 1, r = 0, • ■ • , m — 2, 

lim fl r) (t)^ lim ef (-A) r e A ^b=(-l) r %A r b = 0, 

t^t k - o * t-+t k - o 

by virtue of (2.7). Hence f k (t ) is m - 2 times continuously differentiable, and so is u(t). 

m 

(ii) Set = 0 for t € [0,T], i.e., 

t=i 

m _ 

F(t) = = 0, 

1=1 

m 

with p = ^7 Then we have F^(t) = 0 on [0, T], especially 

i=l 

F {r) (t k ) = p T (-A) r b=0, r = 0, - ■ • m — 1. 

Therefore, 

p(b, Ab, • • • , A m ~ l b) = (F. 

In light of (2.7 ),'p is a zero vector ^ind consequently, 7,- = 0, i = 1, • • • ,m. So <7j’s are m 
linearly independent functions. 

Next we set 

Pn-lfn-l (t) + E7tf«(0 = 0 > te M* (2,U) 

1-1 
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If /?„_! ± 0, then 

• i m 

f n -l(t)=— ^ £ 7 i9i(t)- ( 2 - 12 ) 

Pn-1 i = 1 

m t 

By the definition, f n -i(t) = 0 on So ^7.^(0 = 0 for f € (f n -i,*n) which yields 

1=1 

7, = 0, i = 1, • • • ,m, and hence f n -i(t) = 0 by (2.12). This is a contradiction. Then (2.11) 
yields /?„_ 1 = 0, and consequently 7; = 0, i = 1, • • • , m. So /„_ 1 and gd s are m + 1 linearly 
independent functions. 

Continue the above procedure by adding /jt’s one by one, we are able to show that /V s 
and g^s aure m + n — 1 linearly independent functions. 


3. Construction of splines by the control theory* 

Theorem 2.4 implies that an optimal control for the system (2.1) (with A given by (2.2)) 
is unique. But in general, an d 7 i ? s in (2.10) are difficult to find, we then introduce 
a practical procedure to construct a control law that satisfies all the requirements. This 
control law actually leads us to a construction of spline functions. 

By the existence of a control law, there exists a set of points x 1 ,--*,? 1 1 with x^ = 
a*, k = l,***,n — 1 such that the solution of the system ( 2 . 1 ) satisfies x(£fc) = x fc , k — 
0, 1, * • * , tl — 1 , n. By virtue of Theorem 2.3, a control law that satisfies all the requirement 
can be defined piecewise as 

— u k{t), k = 1 , * • • , n, (3T) 

where ujt(0 is given by (2.9) with t — i*_i, i = fjt, Pl = 1 ? aQ d pR = x* . Then equations 

to find (n — l)(m — 1 ) unknowns in x 1 , • • * , x "” 1 (recall that x* = a* k = 1 , * * ■ , n — 1 , are 
known) come from (n — l)(m — 1 ) continuity conditions on u(i), i.e., 

«*’(**) = «*+,(**). r = 0, ...,m — 2, t=l,...,n-l. (3.2) 

From (2.9), 

„*((») = (A'bfc-^d"' e~ At bS r e~ AT, ds)~ 1 (e~ Atk x k — e~ Atk ~' x k ~ 1 ); (3.3) 

Jt k _\ * 

= («)V ir ‘[ (“*■ - <■'“4 (3-4) 

* tk 

Next, we shall simplify (3.3) and ( 3 . 4 ). Toward this end, we introduce a change of variable 
5 = tk-i + s' into 

( [ tk e-^bFe-^’ds )- 1 = (e - ^*- 1 /** ds'e ~ ATtk - 1 )~ 1 

Jtk~ i . 
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= e ATtk -' M(h k )e Aik -' , (3.5) 

where M(h k ) is defined by (2.8). Substituting (3.5) into (3.3), we have 

u[ T \t k ) = - i*' 1 ). (3.6) 

Similarly, 

u^ih) = (A r b) T M (h k +i)(e~ Ahk+l x k+1 - x k ). (3.7) 

Substituting (3.6) and (3.7) into (3.2) yields a linear system for (n - l)(m - 1) unknowns in 

fV-.,*"- 1 : 

—{A r b) T e~ AT hk M{hk)x k ~ l + ( A T b) T [e~ AThk M{h k )e~ Ahk + M(h k+1 )]x k 
-{A r b) T M{h k+l )e- Ahk +'x k+1 = 0, r = 0, • • • , m — 2, k = 1, • • • ,n - 1. (3.8) 


By virtue of the existence and uniqueness of the optimal control, the linear system (3.8) has 
a unique solution and hence its coefficient matrix is invertible. 

In order to solve (3.8), The following quantities needs to be calculated, A r , e~ Afl (e~ A K ), 
and M(h). Sometimes it is easier to use the Jordan matrix of A, denoted by A. There exists 
an invertible matrix Q such that A = and hence 

A r = QA r Q~\ e~ Ah = Qe~ Ah Q~\ e~ ATh = Q- T e~ ATh Q T . (3.9) 

M(h) = (£ Qe- A 3 Q- l bb T Q~ T e- ATa Q T ds )- 1 = Q~ T M{h)Q~\ (3.10) 

where 

M(h) = {j\-^Q- l b{Q- l b) T t- KTa ds)-\ (3.11) 

Substituting (3.9) and (3.10) into (3.8), we then have 

-(A r Q- x b) T e- KThk M(h k )Q- l x k - 1 + {A T Q- l b) T [e- AThk M(hk)e~ AKk + M(h k+ 1 )}Q~ 1 x k 
— (A r Q~ 1 b) T M(hk+i)e ~ xhk+1 Q~ l x k+l =0, r = 0,---,m-2, k= 1, • • ■ , n - l.(3.12) 


Solving (3.8) or (3.12) for (n - l)(m - 1) unknowns in x l , ■ ■ • ,x n_1 , we then have the control 
it(/) defined piecewise by (3.1). The solution of the system (2.1) is thus given by 


:(0 = 


At -7O , 

e x + 


/‘ e A ^' a) bu{ s ) ds 

Jo 

Q(e M Q~ 1 x°+ f e fi ^ t ~^Q~ 1 bu(s)ds). 
Jo 


(3.13) 
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Note that i'(t) = x,+i(f), * = 1, • • ■ ,m — 1. So the continuity of xj(f) is continuity of 
z, + i(t) for i < m. Further, continuity of Zj m+r ^(<) is continuity of u* r) (t), r = 0, • • • , m — 2. 
Therefore, the observation function y(t) = (Fxft) = Xi(t) is a 2m — 2 times continuously 
differentiable function that satisfies the boundary conditions 

t/ (r) (0) = x° r+11 yU(T) = xj +l , r = 0, • • • , m — 1 (3.14) 

and the interpolation conditions 

y(t k ) = x k l , k = 1, ♦ • • , n — 1. (3.15) 

Hence y(t) is a spline function. We see that from the control theory, we can derive quite 
general spline functions. Summing up, we have proved 

Theorem 3.1 : (1) There exists a unique function y(t ) £ C m ~ 2 [0,T] that satisfies the 
boundary conditions (3. If) and the interpolation conditions (3.15); (2) y(t) is the first com- 
ponent of the vector function x(t) given by (3.13) in which u(s) is defined piecewise on each 
subinterval tk], k = 1, * , n, by 

u{ r \t k ) = l?'e- AT ( t ~ tk ~ 1 ^ M(hk)(e~ Ahk x k — x k ~ l ) 

= (g“ 1 6) r e“ Ar ( t * tfc * l) M(^)(e- A ^Q- 1 f* - Q^x*" 1 ), (3.16) 

where x* , k = 1, • • • ,n — 1 are determined by solving the linear systems (3.8) or (3.12) (Note 
that x°, x 71 and x\, k = 1, • « • , n — 1 are given by the boundary conditions (3. If) and the 
interpolation conditions (3.15)). 

In the next section, we will see that these splines can be piecewise polynomials, trigono- 
metric functions, exponentials or any combination. As special cases, we are able to recover 
classical polynomial splines (odd order) and exponential splines by properly selecting pa- 
rameters ai, * * • , a m in (2.2) for the matrix A. 

4. Classification of splines. 

The type of the splines is determined by its nodal shape functions. From the control 
theory, we are able to construct the nodal shape functions of splines. 

In order to see the kind of interpolation functions in x(f), we only need to consider one 
subinterval. Without loss of generality, we use the first interval (fo?^i) = (0, /t) where the 
solution of the system (2.1) is given by 

X (t) = e At x° + j e A( ' t ~ 3 ^bu(s)ds. (4.1) 
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From Theorem 2.3, 


u(t) = P"e~ ATt ( I* e-^bP'e-^’dsy'ie-* 1 '! 1 - x°). (4.2) 

Substituting (4.2) into (4.1), we have 

x(t) = e At x° + e^'-^MPe-^’dsMihXe-^x' - x°) 

= Qe Xt {Q~ l ? + M(t)- l M{h){e~ Ah Q- l x l - Q~ l Z°)}, (4.3) 

where M(h) and M(h) are defined by (2.8) and (3.11), respectively. 

Theorem 4.1 : Let A be given by (2.2), let (pi(f), — ,PmW) the first row of the matrix 
e At [I - M(t)~ l M(h)] = Qe M [I - M (f)~ l M (h))Q~ x , (4.4) 

and let ( qi(t ), • • • ,q m (t)) be the first row of the matrix 

e At M(t)- l M(h)e~ Ah = ge A ‘M(f)' 1 M(/i)e- AA g- 1 . (4.5) 

Then for r = 0, • • • , m — 1, 

p! r) ( 0) = <W, p! r) W = 0, i = l,---,m, (4.6) 

^ r) (0) = 0, q\ r) {h) = Sj.r+u j = 1," • ,m. (4.7) 


Proof : From Theorem 2.1 and (2.7), the system (2.1) is controllable. By virtue of 
Theorem 3.3, a control that moves x(t) from x(0) = x° to x(h) = x 1 is given by (2.9) (with 
t = 0, t = h, pl = x°, Pr — x 1 ), and consequently 


y(t) = efx(t) 

= %{e M P + f e A(t - s) bu{*)t*) 

= ef e At [I - M{t)- l M(h)}3? + ef e M M{t)- x M{h)e~ Ah x l 

= (pi(0» >Pm(0)®° + (?lW> " ' »9»»W)* 1 

m tn > 

* i=i • j=» i" • 

Choose x° = e,- x 1 = 0 in (4.8), and we have 

p\ r \t) = y (r) (*)> r = 0, "*,m — 1. 


(4.8) 
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Therefore for r = 0, • • • , m — 1, 

p! r >(0) = »<->(0) = x< r, (0) = x.+i(0) = ij +1 = *,«, 

Ppih) = y fr) {h) = l\ r '(h) - x r+ i(A) = x‘ +l = 0, * = 1, 

Then we have proved (4.6). The proof of (4.7) is similar. ■ 


We call pi, qi nodal shape functions by the characteristics (4.6) and (4.7). From (4.4) and 
(4.5), We see that the nodal shape functions are linear combinations of function entries of 
matrices e At and e Ai M(t)~ l . In order to see the type of functions in the spline, we only need 
to examine the entries of these two matrices. 

In the following, we classify the spline functions derived from control theory. This clas- 
sification is based on the spectrum of the coefficient matrix A of the system (2.1) under 
different circumstances. We shall concentrate on the case m = 2. The reasons are: (1) The 
general situation for large m is very complicated and is difficult to describe precisely. (2) 
The case m — 2 has almost all features for the general case. (3) From the practical point of 
view, the case m = 2 is the most useful and important case. Let 


A = 



j3,7 e R 1 . 


The eigenvalues of A are Aj = 7 + -v/7 2 + A 2 = 7 — Vl 2 + P- 


1. -y 2 -j- ft > 0. There are two distinct real eigenvalues. Then the Jordan matrix A of A , 
the transformation matrix Q and its inverse are given by 


A = (o :). «=u A 1 ,)' e- 1 = aT^a : ( X ‘1 1 ) 


Then 


e M = 


e Al< 0 

0 e A2t 


(4.9) 


(4.10) 


From (3.11) (by changing h to t), we have 

^W" 1= (A“A tfio ( % e-°^)( l) (-1,1) ( 


e -Al * 0 

0 e _A3 * 


ds 


1- 


. (1 - e- iAlt )/2 A| ( e -( Al+Aj )‘ - l)/(Ai + A 2 ) 


(A 2 - A,) 2 V V (Al+A2)t " l)/(Ai + A 2 ) (1 - e" 2A ’‘)/2A 2 


(4.11) 


and hence 


1 ( (e A " - e- x '‘)/2A, (.**■' - e A ‘')/(A, + A,) ^ u 

e ~ (A, - A,) 2 ( («-»•' - e A »‘)/(A, + A,) ( e J >' - )' l ' 1 
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l.a. 7 / 0, j3 / 0. In this case Aj, A 2 , — Ai, — A 2 are all distinct. We then have the 
exponential spline with basis functions given by linear combinations of e Al ‘, e -Al ‘, e A2 ‘, e _Alt . 

l.b. 7 = 0, 0 > 0. In this case Ai = — A 2 = the basis functions in l.a. degenerate. 
However, by applying the following limits 


e A 2t — e A,t 


lim 

i Aj + A 2 
e -Aii _ e A a t 

lim — - = lim 

Aj— *> — A* Aj -f- A2 A1+A2— ►O 


. . = lim 

Aj— *— Ai Ai + A2 Aj+Aj-h *-0 

e _Ali — e A,< 


e Alt ( e - ( A 1+ A^ ( -l) _ Xlt 
Ai + A 2 

e _Al< (l — e( A « +Aj )‘) 


= —te 


-A,t 


At _ A | ' “ c ' —te x ' 1 A 

e “ TTT I _^-Ai t ( 0 \ x t _ e - A it)/2Aj j • 


(4.13) 


Ai + A 2 

(4.12) becomes 

_1_ / (e Al * — e -Al ‘)/2A 1 
4A 2 V -te~ Alt (e A 

Hence we have the exponential spline with basis functions given by linear combinations of 
C V& C -V&, i evfc, te-vfr 

l.c. /? = 0, 7 ^ 0. In this case, Ai = 0 (if 7 < 0), or A 2 = 0 (if 7 > 0). Again the basis 
functions in l.a. degenerate. Assume that Ai = 0, then A 2 = A = 27. From the limits 


lim 

Ai— 0 


,Aif Ait 


2Ai 


= t , lim e Alt = 1, 

Ai— 0 


we have 


» At =^ °AtV = Xt 


.-At 


- 1 


0 e 

/ 

e A, M(()-‘ = i 1 X ‘ 


p l - 1 (1 - e-™)/2 ) ’ 


e- At - 1 \ 

» At -e - At )/2 J 


(4.14) 

^ 1 -e A< (e Al - An/ ° 1 • ^ 4 ‘ 15 ^ 

Therefore we end up with the exponential spline with basis functions given by linear com- 
binations of 1, *, e 27 ‘, e~ 2 ' ,t . Later we shall further show that this is the classical exponential 
spline [9]. 

2. 7 2 + j3 < 0. There are two complex eigenvalues: Ai = 7 + iu, A 2 = 7 — z’u;, where 
u = \/— 7 2 - (3. 

2. a. 7 ± 0, fi < 0. Evaluating (4.10) and (4.12), we have the exponential- trigonometric 
spline with basis functions given by linear combinations of e'^sinuif, e 7 ‘coso >t, e ^ sin wt, 
e - ' 1 '* cos cot.,.* * 

2.b. 7 = 0, /? < 0. Again this is a degenerated case where Ai = A 2 = tu> = 
Therefore (4.10) is now 


0 At _ 


coswf + i smut 

0 


0 

cos ut + z smut 


)■ 


(4.16) 
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Taking the limit 7 — ► 0 in (4.12), we then have 

— 1 ( sin ut/uj 


At 1 1 oiuwbyw i.^wuwv 1 1 smart) \ 

e ( / — 4^2 l — t(cosojt — i sinutt) sinwf/u; f J 


— f(cosu>i + 


(4.17) 


Hence, we have the polynomial-trigonometric spline with basis functions given by linear 
combinations of sin cos \/— /?t, t sin t cos 

3. 7 2 + j3 = 0. In this case = A 2 = 7. 

3. a. 7^0. We have non- degenerated Jordan form in this case, 


A = 


7 1 

0 7 


0 - ( \ "1 ) ' (413) 


Therefore, 


e A< = 


e yt te 

0 




)■•*(• 0 - 


(4.19) 




(1/t — ^) 2 1/7- 


- j£ 


i _J ) 




(4.20) 


1/7 -5 

After some more detailed manipulation (see the next section) we can show that this is the 
exponential spline with basis functions given by linear combinations of e 1 , tt 1 , e ^ , tc , 
similar to the case l.b. 

3.b. 7 = 0. In this case, A itself is a Jordan matrix. We compute directly the following 
quantities: 


A = 


0 1 
0 0 


_ At 




M( ‘ r ' =1(5 i 3 )(i) (0,1) ( - ? ) ■ d3 = ( % _! ' /2 

= ( S', l ‘ T ) ■ 


(4.21) 

(4.22) 

(4.23) 


Then we have the polynomial spline with basis functions given by linear combinations of 
* 1 ,t,t 2 ,t 3 . In the next section, we shall further show that this is the well-known cubic spline 

[4]- 

From the above discussion, we see that we may encounter all kinds of splines by varying 
parameters and 7. Two general cases are l.a. and 2.a. where we have full sized exponential 
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or exponential-trigonometric splines. Degeneration occurs when zero or multiple eigenvalues 
appear. The extremal is the case 3.b. when both eigenvalues are zero. It is this extremal 
case that draws most of the attention. This is evidenced by extensive investigation regarding 
the cubic spline in the literature. Case l.c. also has been investigated from a different point 
of view. But we can hardly find any work regarding the other cases (except l.c. and 3.b.) 
listed above. 

The situation for m > 2 is similar. Let Ai, • • * , A m be eigenvalues of A . When Ai, • • • , A m ; 
— Ai, • ■ • , — A m are all distinct, we have the exponential spline with the basis functions given 
by linear combinations of 

e Al *, e~ Al *, •, e Xm \ e“ Amt . 

See case l.a. When complex eigenvalues appear, we get the exponential-trigonometric splines 
with basis functions e Afct sinu;t, e~ Afc£ cosu;i (see case 2. a.). If we have multiple eigenvalues, 
the terms like 

<e A£ , tsinut, te~ Xt cosuit, f 2 e At , 

will appear in basis functions (see cases l.b., 2.b. and 3. a.). Finally, zero eigenvalues will 
introduce polynomials into basis functions (see case l.c.) and the extremal situation is that 
all eigenvalues are zero in which case we recover polynomial splines of order 2m — 1 (see case 
3.b.). 

5. Examples of splines. 

In this section, we shall work out in detail some classes of splines. We shall explicitly 
construct the nodal shape functions and the linear system needed to solve for the unknown 
parameters. 

1. Our first example is the case 3.b. which turns out to be the classical cubic spline. We 
first construct the nodal shape functions. From Theorem 4.1 we need only to calculate the 
first row of matrix (4.4) and the first row of matrix (4.5). Let t = h in (4.22), and we have 

Mth\ — ( *73 -*72 Y'-'K * * 7 *\ 

- ( — *72 * ) ~ h'\ *72 *73 j ' 

Thus from (4.21) and (4.23), we have 
e At — e At M (t)* 1 M (h) = 


M tm*L Hm » 

at ran quftjrr 


t\ ( 1 t \ ( -* 3 /6 *72 \ 12 ( h h>/2 \ 

l )~\0 1 ) { -*72 * ) k* V *72 *73 ) 
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■ (it)" 


—2t 3 /h 3 + 3t 2 /h 2 t(-t 2 /h 2 + 2 t/h) 


)■ 


Hence 


= (-*7* + **/** *€-*»/*« + »/*} J ^ 1 -Aj 

/ —2t 3 /h 3 + 3t 2 /h'‘ l(i 2 /A 2 - 1/A) ^ 

*(() = ! + 2( j) 3 - 3(|) 2 , M ‘) = ‘[1 + (£) 2 - 2(j)l; 

?i(0 — — 2(^-) 3 + 3(— ) 2 , ? 2 (*) = ‘(j)l(^) — 11- 


(5.1) 

(5.2) 


They are precisely the nodal shape functions for the Hermit interpolation. 

Next we set in (3.8), r = 0, x k — (ajt,^jt) T , and h * = Afc+i = A (by this, we are using 
equally spaced intervals). The equation (3.8) is now 

- &e- ATh M{h)x k ~ l + P~[e- ATh M(h)e- Ah + M(h))x k - F M(h)e~ Ah x* +1 = 0, (5.3) 

for k = 1, • • • , n — 1. Substituting 

_ / 1 0 A 12 | 

-h 1 ) A< l A 2 /2 A 3 /3 i 


A A 2 /2 \ _ 12 / 1 A/2 

A3 —A/2 — A 2 /6 


W(A) = 

M{h)e~ Ak = [e-^M(/,r = | ( ,) 2 1$ ) . 


e- ATh M{h)e~ Ah = ^ 


12 


A 3 


1 A/2 

-A/2 -A 2 /6 


W 1 -A \ 12 / 1 -A/2 \ 

J V 0 1 ) A 3 V -A/2 A 2 /3 ) 


into (5.3) yiel’ds, 


6 2.1 a k - 




or 


< SM>VA 


fik-i + 4/3* + A+i = ^(o* +1 “ A = 1, • • • , n - 1. 


(5.4) 


In the matrix form (5.4) is 


/ 4 1 0 ••• 0 ON 


A ^ 


0-2—00 ^ 


A \ 

1 4 1 ••• 0 0 


A 

• 

a 3 -ai 


0 

0 1 4 ••• 0 0 


A 

3 

o 4 — a 2 


0 

: : : : : 


; 

“ A 

l 


• 

0 0 0 ••• 4 1 


fin-7 


On— 1 O n — 3 


0 

f— H 

o 

o 

o 


\ fin-1 J 


On - 0„_2 J 


fin / 


(5.5) 
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This is precisely the same linear system as we construct the cubic spline. After solving /?*’ s 
from (5.5), the desired cubic spline can be expressed piecewisely on the subinterval <*], 
k = 1 , • • • , n as 

y(t) = Q k -lPl(t - tjfc-i) + #fc-lP2(i - tk- 1 ) + oc k qi(t - tk-1 ) + - tk- 1), 

where Pi,P2,9i,92 are given by (5.1) and (5.2). 

2. The second example is the case l.c. which is the classical exponential spline [9]. The 
nodal shape functions axe again derived from the first rows of matrices (4.4) and (4.5). We 
use the Jordan form, let Ai = 0, and A 2 = A (4.9), and we have 

Using (4.14) we can calculate (by setting t — h ), 

= C(A) ( (1 + \ 

2A 3 


M(h) = A e _ A< 


A t e~ Xt - 1 


-1 


a -AA 


where 


-1 (1 — e~ 2Xt )/2 

C( A) = 


)/2 1 

Xh(l - e~ Xh )~ l 


AA(1 + e~ Xh ) — 2(1 — e~ Xfl ) 
Recall (4.14) and (4.15), and we have 


Qt Kt [I - M{t)~ x M{h)\Q~ l 

( 1 1 \( 1 0 )r/_£W 

( 0 A ) [ 0 e AI ) A3 


\ f p~ Xt _ 1 

1 — e Xt (e Xt -e~ Xt )/2 


(l+e~ AA )/2 1 W A -l\ 

.1 Xh(l -e~ xh )~ l ) l X \ 0 1 ) 

1 (e A£ — 1 )/ A \ C( A) 


0 e A£ J 

(1 + e~ Xh )/2 

1 Xh(l 

Pi(Mt) P2(A;t) A 


A 4 


l 1 1 \( A t e~ xt - 1 

^ 0 At ) { 1 — e At (e A< - e- A< )/2 


-e- xh )-' ) ( 0 l 1 ) 


where 


Pi(A;t) = 1 - 


At(l + e~ A/> ) - (1 - e~ A/l ) - e A <«-*> + e~ At 
Xh(l + e~ xh ) — 2(1 — e~ Xh ) 


(5.6) 
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p^(X)t) = 


e A ‘-l 


A 1 + e~ Xh - 2 


l-e —*' 1 l 


1 - e~ Xh t l-e A * 1 - (1-e-Afc) e A * - 2 + e Af 
vT - * n / 


AA 


AA v /i A h 


Xh 


Xh 


]• 


(5.7) 


“)/2 ) 


ge A ‘M(i)- 1 M(/i)e- AA g- 1 
C(A) / 1 1 \ / 1 0 W A t e~ xt 

~\0 aJI^O e xt ){e- xt -l (1 — e -2An 

f (l+e- AA )/2 1 W 1 0 W A -1 \ 

^ 1 AA(1 — e -AA ) -1 ) V 0 e' AA ) ^ 0 1 j 

(7(A) ( 1 1 \f A t t~ Xt — 1 ^(l+e- AA )/2 e- AA 

1 A h(e Xh - l)" 1 


a 4 V 0 A A 1 - eA< ( eXt - e_A ‘)/ 2 

0 v) 


_ ( q\{X-,t) q 2 (X;t) 


where 


qi{X;t) = 1 -pi(A;f), 


(5.8) 


?2(A; t) — 


1 4- e -xk __ 2 i =g^ xh K h Xh 

A/i 


1 p-AA ^ 1 _ pAt A* 1 -At _ 9 , p-Ai 

1_e fl + i — - . lli— ]. (5.9) 

M. 1 \h Xh ' 

(5.6) - (5.9) are nodal shape functions for the classical exponential spline. 

Next we set in (3.12), r — 0, x* = (afc,/?fc) T , and hi t = hk+i = h. The equation (3.12) is 
now 

-(Q- 1 b) T e~ ATh M{h)Q- 1 x*- 1 + (g- 1 6) r [e- ATA M(/i) e - AA + M{h)]Q~ l x k 
—(Q~ l b) T M(h)e~ Ah Q~ l x k+1 = 0, k = 1, • • • ,n - 1. (5.10) 


Substituting 


(Q~ l b) = t(“1’1)’ 


-Ah 


= C i A )(5 e-“)( (1+e i 

Xh(e 

-AA _ r--AA mymiT _ ni \a ^ 


)/2 '* 1 

AA(1 - e _AA ) _1 


- mu ( 1 + e_A '*)/ 2 1 ^ 

- °( A M p-AA \h(fi XX — 1) — 1 /’ 


M{h)e~™ = [e- Ah M{h)Y = (7(A) 


.-aa \ 

AA _ i\-l > 


AA(e AA -l) 
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e~ Ah M(h)e~ Ak = C(A) 


1 0 


(1 + e~ Xk )/2 e~ Xk 

1 A h(e Xk - l)" 1 


- C(X)( ( 1 + e ~ Xh )/ 2 e~ Xk ) 

- G (A)^ C -A* Xhe~ xk {e xk - I)’ 1 ) * 


into (5.10), canceling C(A)/A 2 , we have 

/ ! d ( (1 + e ~^)/2 1 

1 Mj l e" A/l Xh(e xk — l) _l 


A -1 \ ( Qk-\ 

o i M ft-i 




1 + c 


1 + e~ Xk l + e" M W A -1 W ct k 

1 + e~ XK \h{l + e~ 2Xh {\ - e~ Xh )~ l )\0 1 )\ 0 k 

(1 + e~ xk )/2 e~ xh \ / A -l \ ( ojt+i \ 


/ , n Ml + e -**)/2 

■( ’ M 1 Xh(e xk — l) -1 


0 1 J\0k+ 


1 _ e ~ 2AA - 2Xhe~ xk , a Q . .\h(l 4 

2( - r _- e _ A/l) (&-i + &+i) + I— fZ" 

1 - e~ A/l 

= A (ofc+i — a*-i), fc = - 1. 


)]A 


(5.11) 


This is the linear system for the exponential spline, it can be written in the matrix form as 

/ a(X;h) b(X; h) 0 ••• 0 0 W ft \ 

6(A; h) a(A; h) 6(A; h) 0 0 ft 

0 6(A; h) a(A; h) 0 0 ft 


/ a 2 — a 0 
a 3 - ax 

3 ^ 4 - 0:2 


a(A; h) b(X-,h) ft_ 2 
6(A;A) a(X;h) ) \ ft_ a 


6 (A;A)ft 


(5.12) 


0,1—1 O n _ 3 
Oa.— On-2 


4(A;A)& 


where 


.AA(l+e-»»)-(l-e-,») l - - 2Ate~ 

a(X,h)-6 AA( !-«-«)» ’ b ( X ’ h ’- 3 \h(l - e -«) 2 


Again, after finding /Vs, the spline can be expressed piecewise by the nodal shape functions. 
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It is interesting to examine the limiting 
Applying the L’Hospital’s rule three times, 


case for the exponential spline obtained above, 
we are able to verify that 


lima(A; h) = 4, lim 6(A; /i) = 1. 


(5.13) 


We then recover (5.5), the tridiagonal systems for the cubic spline. Further we can verify 
that (by successively using the L’Hospital’s rule) 

limgi(A; h) = - 1 - (^ - l) 3 - (^) 3 = -2(-) 3 + 3(-) 2 = ?i(<)> 


<7x(t) is one of the nodal shape functions for the cubic spline given by (5.2). The other three 
nodal shape functions can be verified similarly. So the cubic spline is the limiting case for 
the exponential spline 4.1.c. when A — * 0. This is not a surprise from the following limit 
regarding matrix A of the system (2.1), 


lim 

A— o 




(5.14) 


where the left hand side is the matrix A for the case l.c., and the right hand side is the 
matrix A for the case 3.b. We can also examine the limit when A — » co where 


lim a(A; k) = 6, lim 6(A; h) = 0, (5.15) 

A— »oo A— *-00 

lim pi(A; t) = 1 — lim <?i(A; i) = — , (5.16) 

A— *00 n A— >oo fl 

lim p2(\;t) = 0 = lim 92(^1 0 * ( 5 - 17 ) 

A— 00 A— *oo 

Substituting (5.15) into (5.12), we then have 0k = (ojt+i — <*jk-i)/2A , k = 1, • • • , n — 1, which 
is the central difference scheme. We see that (5.16) gives the nodal shape functions for the 
linear interpolation. In order to verify (5.17), it is convenient to rewrite 


Pa(A;t) 


g 2 (A;<) 


,-A h 


\h — 2 1 — t~ Xh (\h — 2)[AA(1 + e 


1 C(A), e- M -2 1— 


A + 
C(A) 
2A 3 


2A 3 


[-f(l- e - AA ) + 


(\h + 2)e~ Xh 2 h 

-M) - 2(1 - e~ Xh )y \ 1 — e~ xh 

e~ Xt h 


k )] 


1 


— (e— -2)], 


H( 1 - e~ Xh ) + 


1 _ e -*‘ e _AA - e x (*~ h ) h(e x ^~ h ^ - 2e~ xh + e -A(£+/l) ) 


1 — e 


-A/i 




So the linear spline (the piecewise linear interpolation) is the limiting case for the exponential 

spline 4 . l.c. when A — > 00. 
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3. The third example is the case 2.b. when 0 = — u> 2 , 7 = 0 and 



Denote 

C\ — u 2 h 2 — sin 2 u>h } C 2 = 7: sin 2 u>h + uth, C 3 = sin 2 u;/», C 4 = - sin 2 uih — uh. 

2 * 

Following the same procedure as the first example, we have the nodal shape functions as 
following 


Cipi(t) = Ci cos u>t + C 2 (sinu>t — cut cos ojt) — C3u;t sin ut, 

CiP 2 (t ) = u>h 2 s'mut - C 3 icosud + C^tsinut, 

qi(t ) = 1-Pi(<)» 

C\q 2 {t) = hsmuh(sin.ujt - utcosujt) - tsinu}t(sinuh - uhcosuh). 

When evenly spaced intervals are used, the tridiagonal system for unknowns 0k is given by 
b(u ; h)0k- 1 + a(w; h)Pk 4- b(ui; h)0k+\ = c(w; /i)(orfc+i — ajt-i), k = 1 , • • ■ , n, (5.18) 


where 

a(u;; A) = 2 uh — sin 2 uh, b(u>] h) = sinuh — uihsinuh, c(u?;h)=u 2 hsmuh. 


It is easy to verify that, 


0 [yjh) z 


= 4, 


3 b(uj;h) 
lim , 
w— 0 (u>/i) 3 


= 1 , 


W— 0 


3 

h' 


which are the coefficients for the tridiagonal system of the cubic spline. We can also verify 
that at the limit u — ► 0, the nodal shape functions p,(u;;t) and t), i = 1,2, have the 
relative nodal shape functions of the cubic spline as their limit when u — * 0. Indeed, we 
recover the cubic spline from this trigonometric spline when u — ► 0. 


4. The fourth example is the case 3. a. where 0 — — 7 2 and 



Denote 

C = 1 - 2e- 2 ^ + e- 4 ^ - 4 7 2 /iV 2 A 
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Computing the first rows of matrices e At — e Ai M(t ) 1 M(h) and e At M(t ) 1 M(h)e Ah , we 
have the nodal shape functions 

C Pl (t ) = e r(2A+t) ( 2 -y/i - 2 7 2 /l* - 7 f + 2 7 2 Af - 1) + e" 7 ‘(l + 7 f) 

+e -r(a*-*)( 7t + 2 7 2 / l i _ 2 jh + 2 7 2 h 2 - 1) + e- n,(4 '*- t) (l - 7 t), 

Cp 2 (<) = t{e-' ,{4h ~ t) + e _7 ‘) + e- 7(2/, -‘ ) (2 7 ht - 2 7 h 2 - <) + e -M 2 '>+0( 2 7 / l 2 _ * _ 2 7 /if), 
<7i(0 = l~Pi(0> 

c 92 (f) = e- 7(3A - i >(/i - t - 2-fht) + ( e - 7 < 3A+ ‘> + e- 7(A_t) )(f - h) + e~^ h+t) (2-/ht + h-t ). 


Evaluating (5.3) in the current case, we get the tridiagonal system (5.18) with 

a(o;; h) = a( 7 ; h) = e yh (l — e 47A — A-yhe 2lh ), 
b(u-, h) = 6( 7 ; h) = e~ 2 ~ lh (l + ~fh) — (1 — 'yh), 
c(u>- h) = c( 7 ; h) = 7 2 /i( 1 - e _27A ). 


Again we can verify that the cubic spline is the limiting case for this exponential spline when 
7 — ► 0. 

5. Our last example is a case for m = 3 when 


(0 1 0\ 

A= 0 0 1 . (5.19) 

V 0 0 0 ) 


The system (2.1) with A given by (5.19) produces the quintic spline. In this case, 


e 


At 


M(t)- 1 


M(h) 


( 1 t t 2 / 2 \ 

0 1 t 

10 0 1 / 

/i 5 /( 20) -<7 8 t 3 /6 \ 

f e~ Aa bb r e~ AT 3 ds = -t 4 / 8 t 3 / 3 -t 2 / 2 , 

0 \ t 3 / 6 -t 2 / 2 t ) 


/ 240 /h 5 120 /h 4 20/ h 3 \ 
[M(h)" 1 ]- 1 = 3 120/A 4 64/h 3 12/A 2 . 

\ 20/ A 3 12/A 2 3/k 


( 5 . 20 ) 

♦ 


(5.21) 


(5.22) 


Using (5.20) - (5.22) to compute the first rows of matrices t At — e M M(t) l M(h ) and 
e At M(t)~ 1 M(h)e~ AK , we then have the nodal shape functions for the quintic interpolation: 

Pl (t) = — ~ 5 ^ 3 [h 2 + 3ht + 6£ 2 ], qi (t) = y 5 [h 2 -3ht + 6t 2 }, 
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Pl(t) 

Ps(<) 


+ 3f 2 ], „(<) = £[*(« — h) — 3 ( J ], 

(* ~ <) 3 < ; „ _ < 3 (< - <>f 

2h 3 ' q 3(<) 2/i 3 ‘ 


In order to compute the parameters for the optimal control, we set in ( 3 . 8 ) r — 0, 1 x k — 
(ajt,/?fc,7fc) r > and hjt = hk+\ = h. Then except ( 5 . 3 ), we also have 


- ( Ab) T c- ATh M{h)x k -' + (Abf[e- ATh M(h)e~ Ah + M(h))x k - (Abf M(h)e- Ah x k+ ' = 0, 

(5.23) 

for it = 1, • • • , n — 1. Substituting (5.20) - (5.22) into (5.3) and (5.23), after some tedious 
symbolic manipulation, we have the following linear system, 

20 

8(Pk + 1 - Pk-i) + h(— 7jt_i + 67,- - 7fc+i) = y(/t-i “ 2 A + /*+ 1); (5.24) 


7 fa-i + 16 P k + 7pk+i + A(7 *_i - 7fe+i) = y(/*+i - A-i)» ( 5 - 25 ) 

for it = 1, • • • ,n — 1. If we arrange the unknowns as (/?i, /171, • • • ^7n-i)> we w iM g e t 

a linear system with a banded 6-diagonal coefficient matrix; if we arrange the unknowns 
as (/?!,- ••,/?n-i,/*7 i,"-,* 7n-i) = 0 T ,hi T ), we will have a linear system with a block 
tridiagonal coefficient matrix. 

' 5 E T 
8 E H 

where 


‘ 16 

7 

0 

... 0 

0 ‘ 


' 6 

-1 

0 

0 

0 ■ 

7 

16 

7 

... 0 

0 


-1 

6 

-1 ••• 

0 

0 

0 

7 

16 

... 0 

0 

, H = 

0 

-1 

6 ••• 

0 

0 

0 

0 

0 

• •• 16 

. 7 


0 

0 

0 ••• 

6 

-1 

0 

0 

0 

... 7 

16 _ 


0 

0 

0 ••• 

-1 

6 _ 


f 0 

1 

0 ••• 

0 

0 

r— t 
1 

0 

1 ••• 

0 

0 

0 

-1 

0 ••• 

0 

0 

I 

i 

• ■*. 

l 

j 

0 

0 

0 ••• 

0 

1 

1 

0 

0 

0 ••• 

-1 

O 




/ 

/l 7 


9 J 


(5.26) 
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’ 15(/ 2 - fo)/h - 7fi 0 ~ h~f 0 


20(/o — 2/i + + 8/?o + A70 

15(/3 — fi)/h 


20(/i — 2/2 + fz)lh 

15 (/<-/ 2 )M 

, 9 = 

20(/ 2 — 2/3 + /i)/^ 

15(/n-l - fn-3)/h 


20(/„_3 — 2/„_2 + f n -i)/h 

15(/n — fn- 2)//l — 7 (3 n + h'fn _ 


. 20(/ n _2 — 2/„_i + fn)/h ~ 8 / 3 „ + h~fn . 


In general, if A is a m x m nilpotent matrix with l’s on the super diagonal and 0’s 
elsewhere, we shall recover all odd degree polynomial splines (with degree 2m — 1). 

6. Numerical experiments. 

In this section, we test the behaviors of different splines numerically. Equally spaced 
intervals axe used for all computations. 

Example 1. Comparison of the cubic spline with the quintic spline. 

Test function 1. 

' 0 -1 < t < 0 

f(t) = <j 1/2 t = 0 

1 0 < t < 1 

For the cubic spline, we pose, in (5.5), the boundary conditions: 

fio = 0 = fi n ', 

and for the quintic spline, we pose, in (5.25), the boundary condition: 

fio = 0 = /?„, 7o = 0 = 7n- 

Recall that /?, and 7 y are the coefficients for the first and the second derivatives, respectively. 
The spline functions axe then constructed for k = .2, h = .1, h = .05, and h = .025. Graphs 
are plotted in Figure 1(a), 1(b). We see that the qualitative behavior of the two splines are 
• * almost same, but the quintic spline has a little better accuracy. 

One interesting phenomenon is that the mesh refinement does not effect the maximum 
overshoot of the spline approximation. Since this is very similar to the Gibbs phenomenon 
for the Fourier series, we term it as “Gibbs phenomenon” of splines. In fact, all spline 
functions have this property. 

Test function 2. 

g(t) = e- l0t \ — 1 < * < 0 . 
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For the cubic spline, we pose the boundary conditions: 


Ao = — 30e 10 , /?„ = 0; 

and for the quintic spline, we pose the boundary conditions: 

Ao = — 30e 10 , An = 0, 7 o = 960e 10 , 7n = 0. 

We use the mesh size h = .2, and plot the graphs in Figure 1(c), 1(d). We see that the 
quintic spline gives much better approximation in the neighborhood of x = —1 since it has 
the correct concavity information at x = — 1 which the cubic spline does not have. 

Example 2. Properties of the classical exponential spline case l.c. 

The test function is the same f(t) as in Example 1. We have observed that for small 
parameter A, the behavior is much like the cubic spline. This is not surprise from (5.14). The 
interesting fact is: For the moderate A, the graph is very much the same as the cubic spline 
(Figure 2(a)). If we fix the parameter A and refine the mesh, we observe Gibbs phenomenon 
as in the cubic and the quintic splines. But if we fix the mesh (here we choose h = .1) and 
increase the parameter A, we see that the approximation converges to the piecewise linear 
function (Figure 2(b), 2(c), 2(d)). This confirms our theoretical analysis made in Section 5. 

Example 3. Properties of the exponential spline case 3.a. 

We use the same test function f(t) as in the examples 1, 2. 

For small parameter 7 , the approximating feature of this spline is also like the cubic spline 
including the Gibbs phenomenon. But when we fix the mesh (here h = .1) and increase the 
parameter 7 , an unexpected wiggling appears at t = 0 (Figure 3(a)-3(d)). 

Example 4. Properties of the exponential spline case l.a. 

Here we choose 

A = 

and we have A! = —1, A 2 = —2. Again, the testing function is f(t) as in the previous 
examples. 

We plot the approximation for h = .1, h = .05, h = .025 in Figure 4(a), 4(b), 4(c), 
respectively, again, we observe the similar behavior as that of the cubic spline. 

Conclusions 

1. Gibbs phenomenon exists for all splines. 

2. The quintic spline is recommended if the concavity is important. 
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3. From the approximation point of view, the classical exponential spline is preferred 
when the function has points of discontinuity. 

A final remark. 

For the discussion purpose, we constructed spline approximation in this paper by intro- 
ducing the nodal shape functions which is not necessary in practical computation. From the 
framework we have established based on the control theory in Section 3, all we need to do 
is: providing the matrix A, the vector b to the linear system (3.8), solving (3.8) numerically 
to obtain x*’s, and hence the control law u(f) (see (2.9)). After we have the control u(i), the 
expected spline function is given by the first component of x(£) defined by (3.13). Based on 
our analysis, we are able to choose different splines by simply selecting entries of the matrix 
A. 

The significance of this investigation is two fold: first, it exposes the relationship between 
two important fields - control theory and spline approximations. This enables us to discover 
new spline functions and to investigate, systematically, the properties of the spline approx- 
imations. Secondly, it provides a practical way to construct different splines from a same 
simple framework. From our experience, we feel that this construction is more natural and 
easier than the traditional approach. 
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